SOLVING ELLIPTIC FINITE ELEMENT SYSTEMS IN 
NEAR-LINEAR TIME WITH SUPPORT PRECONDITIONERS* 

ERIK G. BOMANt, BRUCE HENDRICKSON*, AND STEPHEN VAVASIS§ 

Abstract. We consider linear systems arising from the use of the finite element method for solv- 
ing scalar linear elliptic problems. Our main result is that these linear systems, which are symmetric 
and positive semidefinite, are well approximated by symmetric diagonally dominant matrices. Our 
framework for defining matrix approximation is support theory. Significant graph theoretic work 
has already been developed in the support framework for preconditioners in the diagonally dominant 
case, and in particular it is known that such systems can be solved with iterative methods in nearly 
linear time. Thus, our approximation result implies that these graph theoretic techniques can also 
solve a class of finite element problems in nearly linear time. We show that the support number 
bounds, which control the number of iterations in the preconditioned iterative solver, depend on 
mesh quality measures but not on the problem size or shape of the domain. 

1. Introduction. Finite element discretizations of elliptic partial differential 
equations (PDEs) give rise to large sparse linear systems of equations. A topic of 
great interest is preconditioners for iterative solution of such systems. We focus 
on scalar boundary value problems of the form V • (f?Vit) = — /, in which 6 is a 
scalar conductivity field. See (3.1) below for a more detailed statement of the PDE 
under consideration. Such PDEs arise in a variety of physical applications listed in 
Section 3. We prove that the stiffness matrix K of this PDE, which is defined precisely 
by (3.7) below, can be well approximated by a diagonally dominant matrix K . 'Well 
approximated" means that n(K, K) has an upper bound in terms of mesh quality 
measures only; the notation k(-, •) is defined below. Since significant theory has been 
developed for diagonally dominant matrices, our result shows that the same theory 
extends to this class of stiffness matrices. In particular, our approximation result 
means that the system Kx = f arising from the finite element method (FEM) can be 
solved in nearly linear time by preconditioned iterative methods. 

Our analysis uses the support theory framework described in [5] for analyzing 
condition numbers and (generalized) eigenvalues for preconditioned systems. Our 
analysis is as follows. In Sections 5-7, we state and prove our theorem that K may 
be factored as K = A T D 1 / 2 HD 1 / 2 A. A preliminary general-purpose result in Sec- 
tion 2 shows that a matrix of this form can be approximated by K — A T DA (a 
diagonally dominant matrix), with the quality of approximation depending on k(H). 
The analysis of k(H) in Section 7 establishes that this quantity depends on the space 
dimension d and degree p of the FEM, the quadrature rule, and various mesh quality 
measures but does not depend on the number of nodes or elements or on the shape 
of the domain or size of elements. It also does not depend on the conductivity field 8 
under certain assumptions to be made below. 
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The idea of approximating FEM systems by diagonally dominant matrices is not 
new; see for example Gustafsson [13]. In fact, our approach is similar to Gustafsson's 
in that we approximate each element matrix by a diagonally dominant matrix. In 
contrast to [13], we are able to rigorously prove bounds on the spectral properties of 
our approximation (and thus also the preconditioner) . 

An approach in [3, §7.3] proposes to precondition general non-diagonally dominant 
FEM problems with a diagonally dominant preconditioner obtained by using a lower 
order method. Again, our result is a departure from theirs because we get precise and 
rigorous bounds on the quality of the approximation. 

In Section 10, we specialize the theory to four common cases of the finite ele- 
ment method and report on computational tests for these cases. The computational 
tests are intended to illustrate the size of the constants in our main theorem and do 
not measure the practical performance of our solution strategy versus competitors. 
See further remarks on the practical efficiency of our preconditioning approach in 
Section 9. 

2. Support Theory. Let A, B be symmetric positive semidefinite (SPSD) ma- 
trices of the same size, and let N{A) and N(B) denote their null spaces. We define 
the support number of A with respect to B to be 

a(A,B)= sup -5—. (2.1) 

xeR"-JV(B) x BX 

Note that this quantity is finite only if N(B) C N(A). We follow the convention 
established in the previous literature of using a(-, •) to denote support numbers. Un- 
fortunately, a is also commonly used to denote singular values. In this paper, a with 
one argument is a singular value, and a with two arguments is a support number. 

The results in this section can be partially generalized to indefinite symmetric 
matrices. This generalization requires a more general definition of support number 
than (2.1) and requires the use of the Symmetric Product Support Theorem from [5]. 
Since this paper focuses on the positive semidefinite case, we omit this generalization. 

If A and B are symmetric positive definite (SPD), then <r(A, B) = X max (A, B), the 
largest generalized eigenvalue. When B is a preconditioner for A in the preconditioned 
conjugate gradient iterative method [10], the condition number of the preconditioned 
system is given by cr(A, B)a(B, A), which we denote as k(A, B). When A, B are SPD, 
then k(A, B) = k(B^ 1 A) where n denotes the standard spectral condition number. 

A principal goal of this paper is to propose the construction of a symmetric 
diagonally dominant matrix K that approximates the finite element stiffness matrix K 
in the sense that n(K, K) is not too large. The following two lemmas and subsequent 
theorem define our framework for defining K and bounding k(K, K). 

Lemma 2.1. Suppose V e R nxp , and suppose H e W x p is SPD. Then 

a(VHV T ,VV T ) < A max (ff), 

where A max denotes the largest eigenvalue. 
Proof. 

a(VHV T ,VV T )= sup X 7 T ^ TX 

Y T Hy 

= SUp = — 

y eR(v T ) y y 

< A max (iJ). 
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The second line follows from the first by substituting y = V T x and the third from the 
Courant-Fischer minimax theorem [10, §8.1.1]. Here, R{V T ) denotes the range-space 
of V T . □ 

Lemma 2.2. Suppose V G M™ xp , and suppose H e W xp is SPD. Then 
a(VV T ,VHV T ) < l/X min (H), 

where A m j n denotes the smallest eigenvalue. 

Proof. Using the same idea in the previous proof leads to a suprcmum over the 
Rayleigh quotient y T y/(y T -ffy), which is bounded above by l/A m i n (-ff). □ 

Combining these two lemmas yields a result for the condition number. 

Theorem 2.3. Suppose V e R nxp 7 and suppose H e W xp is SPD. Then 

k(VHV t ,VV t ) < k(H). 

Proof. The result follows from Lemma 2.1 and 2.2 and by using the fact n(H) = 
A max (^)/A min (^) when H is SPD. □ 

The way we will apply this theorem is to let V = A T D 1 ^ 2 , where A is the 
node-arc incidence matrix of a graph and D is a diagonal weight matrix. Then 
K = VV T = A T DA is diagonally dominant while K = VHV T = A T D 1 / 2 HD 1 / 2 A 
is not (in general). Therefore, we now have a tool to approximate non-diagonally- 
dominant matrices. 

3. Finite element analysis. In this section we provide a brief summary of 
isoparametric finite element approximation as well as an introduction to notation 
that is crucial for our main theorem. The material in this section is standard [15] in 
textbooks, except that our description herein uses notation for indexing that is more 
detailed than usual. This section concludes with a summary of our notation. 

The class of problems under consideration consists of finite-element discretizations 
of the following second-order elliptic boundary value problem. Find u : £1 — > M. 
satisfying 

V • (0Vu) = -/ on Q, 

u = uq on r l7 (3.1) 
9du/dn — g on r 2 . 

Here, f2 is a bounded open subset of M. d (typically d = 2 or d — 3), Ti and T 2 form 
a partition of <9f2, 9 is a given scalar field on O that is positive-valued everywhere 
and is sometimes called the conductivity, f : il — > M is a given function called the 
forcing function, uq is a given function called the Dirichlet boundary condition and g 
is another given function called the Neumann boundary condition. 

This problem has applications to many problems in mathematical physics. For 
example, u can represent voltage in a conducting medium, in which case 9 stands for 
electrical conductivity. The two types of boundary conditions stand for, respectively, 
a boundary point held at fixed voltage or a boundary point electrically insulated (or 
with prespecified nonzero current). Another application is thermodynamics, in which 
u stands for temperature of the body, 9 for thermal conductivity, Ti for a boundary 
held at fixed temperature, and T 2 for a boundary insulated (or with prespecified heat 
flow). Problem (3.1) is also used to model membrane deflection and gravity. It arises 
as a subproblem in fluid flow modeling. 

3 



For clarity, we keep track of our assumptions by explicitly numbering them. One 
assumption has already been made: 

Assumption 1. For all iefl, 6{x) > 0. 

Without this assumption, the problem may be ill posed. 

The first step in the isoparametric finite clement method is to produce a mesh of 
this domain. For the remainder of the paper, we assume isoparametric elements are 
defined with respect to the usual polynomial basis on a simplicial reference element, 
although the results can be generalized to other reference elements and basis families. 

In more detail, let To denote the standard unit <i-simplex: for d — 2, this sim- 
plex is the triangle with vertices (0,0), (1,0), (0,1), and for d = 3 this simplex is 
the tetrahedron with vertices (0,0,0), (1,0,0), (0,1,0), (0,0,1). Let p denote the 
polynomial order of the finite element method. In the reference element, position 
an evenly spaced array of I = (p + l)(p + 2)/2 reference nodes when d = 2 or of 
/ = (p + l)(p + 2){p + 3)/6 nodes when d = 3. These nodes have coordinates of the 
form (i/p,j/p) (in two dimensions) or (i/p,j/p, k/p) (in three dimensions) where i, j 
(and, in three dimensions, k) are nonnegative integers whose sum is at most p. Let 
these nodes be enumerated zi,...,zi. 

One generates a mesh of Q composed of m elements to be defined via mapping 
functions. For each element t — 1, . . . , m, there is a mapping function (fit that maps 
T to M. d . Let V(/>t denote the derivative (Jacobian) of (fi t - The following assumption 
is nearly universal in the literature. 

Assumption 2. Mapping function <j>t ■ T — » R d is injective, anddet(V(fi t {z)) > 
for all z 6 T , t = 1, . . . , m. 

The tth element is defined to be <fit(To) and hence has a shape of a curved 
simplex. The function <fi t carries the I reference nodes z\ , . . . , z; to I real-space 
nodes (fit(zi), . . . , (fit(zi), which are often just called nodes. These may be denoted 
Ct,i; • ■ ■ Xt,i- The nodes are chosen so that the nodes on the boundary of the tth 
element coincide with the corresponding nodes on the boundaries of its neighbors. 

We restrict (fit to be a polynomial of degree p, in which case it is uniquely deter- 
mined by the positions of its real-space nodes. In more detail, let N^, [i = 1, . . . , I, 
be a real-valued degree-p polynomial function on T with the property that 



It is not hard to write down an explicit formula for iV^ and to prove that (3.2) 
uniquely determines among degree-p polynomials. These functions Ni, . . . , Ni are 
called shape functions. For each t, (fit necessarily has the formula (fit = (t,i-Wi + ' ' ' + 



There are many duplicate entries in the list • • • > Cm,l because of common 
nodes at the boundaries of the elements. Let wi, . . . , w n i be a listing of the real-space 
nodes with all duplicates removed. This renumbering is specified by an index mapping 
function, the local-to- global numbering map, denoted LG and carrying an index pair 
(t, /i) to an index i £ 1, . . . , n' so that u>i = Q tfJi . 

Finally, the mesh T is specified by listing the nodes u>i, . . . , w n i and the local- 
global mapping defined by LG. From this data, one can deduce all of the 0t's. 

Let £1 be the union of the elements, which may be slightly different from f2. This 
is because at the boundary of f2, the elements may either protrude a bit outside ft or 
may fail to cover a small part of the boundary. Let the boundary of f2 be partitioned 
into Ti and T 2 in correspondence with the partition ri,r 2 of the boundary of O. It 




(3.2) 



Ct,iN t . 
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is necessary to transfer 8 (and several other fields) from f2 to f2 and to transfer the 
boundary conditions to r\ and F 2 . We omit the details. 

The next step in the isoparametric finite element method is to define basis func- 
tions 7Ti , . . . , 7r„< , which are functions from f2 to R satisfying 

7T 4 (w fc ) = j J jf.^J? Jfor fc = l,...,n'. (3.3) 

We require that Wi, when restricted to an element t, must be of the form o 
4>t 1 where LG(t,^i) = i or else must be identically if i docs not occur among 
LG(t, 1 ), . . . , LG(t, I). This requirement uniquely determines the 7Tj's. It can be shown 
that TTi is continuous but fails to be differentiate at inter-element boundaries. 

Let the global numbering be chosen so that w±, . . . , w n , the first n real-space 
nodes, are the nodes with no Dirichlet boundary condition imposed, and let the 
final n' — n nodes have Dirichlet boundary conditions. Now we can define the exact 
assembled stiffness matrix K exact of the finite element method to be the n x n matrix 
whose (i, j) entry is given by 

K c ™ ct {i,j) = [ Vtt 4 (x) • e(x)V-JTj{x) dx (3.4) 
Jh 

where V as usual denotes the gradient. 

Matrix K cxact is sparse, symmetric and positive semidefinitc. Symmetry is ob- 
vious; semidefinitcness follows from a fairly straightforward argument that we omit, 
and sparsity follows because if exact (i, j) is nonzero only if there is an element t that 
contains both nodes Wi and Wj . 

Integral (3.4) is difficult to compute directly because evaluating iti requires eval- 
uation of (frt . Fortunately, this difficulty is avoided by breaking the integral into a 
sum over elements then carrying out the integral over the reference domain following 
a change of variables as follows. 

m ,. 

#exact ( . )j)= ^ / V;c 7r 4 (^(z)).^ t ^))V^(^(z))det(V0 t ^))^, (3.5) 

t=l J T <> 

where the notation V x means derivative with respect to the coordinates of element t 
(as opposed to derivative with respect to z, the coordinates of To). 

The integrand of (3.5) is evaluated using the chain rule for derivatives. Assume 
Wi is a node of element t (else the above integral is 0) . Let \i be the index such that 
LG(t, fi) — i, so that Hi = o cj)^ 1 on element t. Then 

= V</> t (z)- T - VJV M (z), (3.6) 

and similarly for iTj . Here, V0t {z)~ T denotes the transposed inverse of the dxd matrix 
V0i(z), which exists by Assumption 2. The '•' notation in the previous formula 
indicates matrix-vector multiplication. Gradients are regarded as length-<i column 
vectors. 

We assume that the entries of K are not the exact value of the integral (3.5) 
but are obtained by a quadrature rule that we now discuss. Let n, . . . ,r q be points 
in the interior of the reference element To called the Gauss points. (As is common 
practice, we use this terminology even if the quadrature rule is not derived from 
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Gaussian quadrature.) Let u>i,...,u) q be corresponding Gauss weights. We denote 
this quadrature rule, i.e., the set of ordered pairs (ri,Wi), . . . , (r q ,cu q ), by the symbol 
Q. Then in place of (3.5) we take 

m q 

K ( l >j) =££v x 7r i (^ t (r fc )) ■6(<j> t (r k ))V x ir J (<f> t (r k ))dct(V<t> t (r k ))Lu k , (3.7) 
t=i fe=i 

in which V x 7Tj(</>t(rfc)) is evaluated by substituting z = r k into the right-hand side of 
(3.6) and similarly for \7 x ^j(4>t(rk))- Symmetry and sparsity of K follow for the same 
reason as for X cxact ; positive semidefiniteness will follow from results in the next two 
sections (under some additional assumptions to be made). 

We close this section with a summary of the notation introduced thus far. Integers 
that define the size of the computation are: 

d = space dimension of (3.1), 

p = polynomial order of the finite element method, 

(p+l)(p + 2)/2 if d = 2, 



I = number of reference nodes 



(p + l)(p + 2)(p + 3)/6 if d = 3, 



to = number of elements, 

n = number of non-Dirichlct real-space nodes, 
n' = total number of real-space nodes, 

q = number of Gauss points in the quadrature rule Q. 

Sets of points in R d include 

Zi, . . . , zi = reference nodes, 
Ci.ii ■ ■ • j Cm,/ — real-space nodes (local numbering), 
wi , . . . , w n ' — real-space nodes (global numbering) . 

The quadrature rule Q is defined by 

ri, . . . , r q = Gauss points, 
ui, . . . ,cu q = Gauss weights. 

Important domains are: 

To — the reference element, 
O = the domain of (3.1), 
0i (T ), . . . , <f> m (T ) = the elements, 

Q = the approximation to ft given by (f>i(T ) U • • • U (f> m (T ). 

Important functions are: 

<pi, ... , 4> rn = clement mapping functions (To —* 
N\, . . . , Ni — shape functions (T — ► K), 
9 = conductivity (Q — > R), 
LG = local-to-global index mapping {1, . . . , to} x {1, . . . , 1} —> {1, . . . , n'}, 
7Ti , . . . , n n i — basis functions (O — > M) . 
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Finally, scalar quantities to be introduced in the next section, but which are 
included here for completeness, include 



m Qi Mq — min and max quadrature weights, 



9{T) = intra-element conductivity variation, 
Qi, . . . , a m = maximum compression in elements, 
Pi, ... , p m — maximum stretch in elements, 



«i(T), k 2 (T) = mesh quality measures, 

<JQ,p,TQ. p — max and min singular values of Sq. 



This section presented in detail the construction of K, the finite clement stiffness 
matrix, but omitted the construction of f , the right-hand side of the linear system 
Kx = f under consideration. The construction of f is not relevant to the main result 
of this paper, but for the sake of completeness, we give a very brief overview and refer 
the reader to [15] for the details. The steps involved in obtaining a formula for the zth 
entry f are as follows. Multiply the PDE on both sides by 7Tj, integrate over O, and 
apply integration by parts. Applying integration by parts causes a boundary term 
involving Neumann data g to enter the equation. Write u as a linear combination of 
7Tj's for j = 1, . . . , n'. The coefficients of this linear combination are the entries of x for 
j = 1, . . . , n and are known from the Dirichlct boundary condition for j = n+1, . . . , n'. 
Finally, gather terms not involving entries of x together on the right-hand side and 
apply quadrature rules to obtain the ith entry of f . Once the linear system Kx = f is 
solved, the approximation of the unknown function u of (3.1) is recovered as a linear 
combination of the nj's using entries of x and the Dirichlet conditions as coefficients. 

4. Condition numbers and assumptions. Our main theorem about matrix 
approximation, which is Theorem 5.1 below, depends on several scalars associated 
with the finite element method and the problem at hand that we define in this section. 
In addition, this section states further assumptions about the problem and method. 

Two constants appearing in our bound are vtlq and Mq, which we define to be 
the minimum and maximum weights in the quadrature rule, i.e., 



Assumption 3. The quadrature weights are positive, i.e., tuq > 0. 

There is some loss of generality with this assumption because a few popular finite 
element quadrature schemes (but certainly not all) use negative weights [8]. 

Assumption 4. The quadrature scheme is exact for polynomials of degree up to 
2p — 2, i.e., if ip : Tq — > K is a polynomial of degree 2p — 2 or less in the d coordinates 
of Tq, then 



This assumption is quite reasonable since it is usually required anyway for accurate 
solution by finite element analysis: one wants accurate quadrature of • V7Tj. 
We now define a dq x (I — 1) matrix Sq, p according to the following formula: 



vhq = min(^i, . . . uo q ); Mq = max(wi, . . . , u q ). 



(4.1) 




k=l 



( VJV 2 (n) 



VN 3 (n) 



VJVi(ri) \ 




(4.2) 



V ViV 2 (r g ) 



VN 3 (r q ) 



VJVi(r,) / 
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Although we introduce this matrix in this section in order to define two associated 
scalars, the motivation for this definition will be postponed until Section 6. Before 
defining these scalars, we require the following lemma. 

Lemma 4.1. Under Assumption 4, matrix Sq iP has full column rank. 

Proof. Let v <E R' -1 be chosen so that Sq. p v — 0. Define the function U : Tq — > M 
according to the formula U — viN 2 + v 2 N 3 + ■ ■ - + vi-\Ni. Observe that, by definition 
of Sq jP , the first d entries of Sq >p v are precisely VU(r\), and the next d entries 
are WU(r2), etc. Therefore, VJ7 vanishes identically at r\, . . . ,r q . Let ip : Tq — > M. 
be defined by ip = VU ■ VU. Then ip has degree at most 2p — 2, is a nonnegative 
function, and also vanishes identically at n, . . . ,r q . By Assumption 4, this means that 
J T ip = 0. But since ip is nonnegative- valued, we conclude that ip = 0. Therefore, 

= also, and thus U must be a constant function on Tq. But the definition of 
U omits the Ni term, and therefore U{z{) = 0. Since U is constant, this means that 
U = on all of To, and in particular, U(z2) = ■ ■ ■ = U(zi) = 0. But U(z2) = Vi, 
U(z 3 ) = V2, ■ ■ ■ , U{zi) — vi-i by construction of U. Therefore, the entries of v are 
all zeros. Since v was an arbitrary vector in N(Sq, p ), this argument proves that this 
nullspacc contains only the vector, hence Sq, p has full column rank. □ 

The next constants that appear in our theorem are o~q, p and TQ,p, defined as 
follows: 



a Q,P — cr max(S'Q,p); TQ iP = a m i n (SQ iP ). (4.3) 

These constants depend only on p and the quadrature scheme. It follows from the 
lemma that both are positive. 

The remaining scalars and assumptions in this section pertain to the mesh. Define 
for t = 1, . . . , m, 

a t - max{||V0 t (n)- 1 ||, . . . , IIV-M^)" 1 !!}, (4.4) 
A = max{||V0 t (n)||,... ) ||V0 t (r,)||}. (4.5) 

The norms in these equations are matrix 2-norms: |ji?|| = cr max (i?). Finally, for the 
whole mesh, a quality measure is 

Ki(T) = max a t (3t- (4.6) 

t — l,...,m 

Although the a's and /3's are subscripted only by t, it is clear from the definition 
that they also depend on Q. On the other hand, it is possible to get a Q-indcpendcnt 
definition of these by simply taking the upper bounds similar to (4.4), (4.5) over all 
z G T instead of just n, . . . ,r q . Unfortunately, for higher order elements (p > 1), 
there is no simple method to compute max z6 T ||V0t(,z) -1 ||; a technique for obtaining 
an upper bound appears in [21]. 

It should be noted that Ki{T) > 1 since for any z, | V^z) -1 1| • ||V0 t (z)|| > 1. If 
all the elements are well-shaped, i.e., not too distorted when compared to the reference 
element, then n\{T) will not be much larger than 1. 

The second mesh quality measure is 

max fc= i ) ... i9 det(V0 t (r fc )) . . 

k 2 T = max — : ; — — — . 4.7) 

t=i,...,m mmti,...,, det(V(Mr fe )) 

This quantity measures the maximum over elements of the variation in volumetric 
distortion over the element. This may be regarded as a measure of how much elements 



depart from linearity (flatness). Measure k 2 is not completely independent from k,\ 
as the following argument shows. It follows from Hadamard's inequality that 

a- d <det(VMrk)) < 

and therefore 

max fe= i i ... ;9 det(V</> t (r fe )) ^ ^ Q ^ d 



min fe= i ) ... ) gdet(V</>t(r fe )) 



< (a t p t ) 



hence ^(T) < Ki(T) d . This bound is not likely to be tight in practice (see our 
computation results in Table 10.3), so we prefer to distinguish the roles of K\ and k 2 . 

The final two assumptions are qualitative in nature and are meant to indicate 
cases in which our method is expected to work well. Violation of these assumptions 
does not invalidate our theorem but merely indicates that our bound on k(K, K) may 
grow large. 

Assumption 5. Mesh quality measures ki(T),K2(T) are not too large. 

In the case of linear elements (p = 1), Ki(T) is proportional to the reciprocal of 
the minimum angle (in two dimensions) or solid angle (in three dimensions) of the 
mesh and K2(T) = 1. 

Assumption 5 does not imply that the elements are of a uniform size: a uniform 
rescaling of element t does not affect the product a t /3 t . 

The final constant and assumption pertain to the conductivity field. Define 

9(1) = max — : — — — — . (4.8) 

t=i,...,m minfc = i ) ... )9 0(</>t(nfe)J 

In other words, 9(T) measures the maximum intra-element variation of the conduc- 
tivity field 9. 

Assumption 6. Intra-element conductivity variation 9(T) is not too large. 

This assumption implies that if there are huge jumps in 9 in the domain (e.g., 
because one is modeling a domain composed of two materials with vastly different 
conductivities), then the mesh boundaries should be aligned with the conductivity 
jumps. If the mesh boundaries are aligned with the conductivity jumps, then no 
element will have large variation in 9 among its Gauss points and thus 9(T) will not 
be large. 

Note that, as in the case of and n 2 (T) 1 scalar 9(T) also depends on the 

quadrature rule, but this dependence can be eliminated by overestimating 9(T) as 

max zeTo 9(4> t {z)) 

max — n/j. i w ■ 

t=i,...,m mm z6To 0{<Pt(z)) 

5. The matrix approximation. Our main matrix factorization result is sum- 
marized by the following theorem whose proof is explained in upcoming sections. 

Theorem 5.1. Let K be defined by (3.7) above, and let Assumptions 1-4 hold. 
Then K may be factored as A T J T DJ A, where 

• A is an (l — l)m x n reduced node-arc incidence matrix of a certain multigraph, 

• J is a dqm x (I — l)m matrix that is well conditioned in the sense that 

(5.1) 

and 

o- m in{J) > tq, p /ki(T). (5.2) 
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• D is a dqm x dqm positive definite diagonal matrix. 

Further, J T DJ may be refactored as D^^HD 1 ! 2 where D is a (I — l)m x (I — l)m 
SPD diagonal matrix and H is a (I — l)m x (Z — l)m SPD matrix whose condition 
number is bounded as follows: 

M Q (T 2 „ 

k(H) < 9{T) Ki (T) 2 k 2 (T) ■ - f p . (5.3) 

m Q T Q tP 

This theorem can now be combined with Theorem 2.3 to obtain a good approx- 
imation K to the stiffness matrix K . In particular, we take K = A T DA. Note that 
this matrix, being a weighted graph laplacian, is symmetric and diagonally dominant. 
Then in the context of Theorem 2.3, V — A T D X / 2 . The theorem thus implies that 
the condition number of K with respect to K depends only on the condition number 
of H for which we have a good bound (5.3). Then K can be preconditioned using 
techniques in the previous literature. The complete description of the algorithm is 
given below in Section 8. 

The approximation bound clearly depends on the quadrature rule, mesh quality, 
and intra-element variation of 9. It is also informative to make a list of quantities on 
which k(K, K) does not depend: 

• The number of nodes or elements, 

• The size of the elements, 

• Variation in the size of elements (i.e., gradation of the mesh), 

• The shape of 0, 

• The conductivity field 9 (provided that the mesh respects internal boundaries 
where large conductivity jumps occur). 

Although the shape of £1 does not directly affect k(K,K), naturally, it affects the 
mesh generation procedure which in turn affects mesh quality. A particular case that 
requires further consideration is when has a sharp corner; a sharp corner of ft forces 
the mesh to contain at least one element with a sharp corner, thus causing K\(T) to 
be large. Even in this case, however, it may be possible to obtain an upper bound on 
k(K, K) independent of this sharp angle using the analysis of Phillips and Miller (see 
Section 9). 

6. The first factorization. In this section we develop the first factorization of 
K stated in Theorem 5.1. This factorization is related to one proposed in [22]. We 
state the main result of this section as a lemma. 

Lemma 6.1. Matrix K defined by (3.7) can be factored as K = A T S T R T DRSA 
where the factors are as follow. 

• Define matrix A £ R(i-i)"»xn ^ ^ e a S p arse ma trix all of whose entries are 
— 1, or 1. In more detail, A is written in block form 

t M \ 

A = \ 

V A m J 

where At is (I-l)xn for each element t — 1, . . . , m. The columns are indexed 
1 , . . . , n in correspondence with nodes w\ , . . . , w n . Row /x — 1 of A t has a I'm 
column LG(t, /x) for /j, = 2, . . . , I and a '—1 ' in column LG(t, 1). Thus, most 
rows of A have exactly two nonzero entries. If LG(t, /i) > n (i.e., node Qt,n 
lies in the Dirichlet boundary Ti), then the 1 ' entry is omitted. Similarly, if 
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LG(t, 1) > n, then the — 1 ' entry is omitted. For this reason, a few rows of 
A have just one nonzero entry or none at all. 
• Define S € Rdqmx(l-l)m by 



Sq, p 



S Q,p J , 



> m times, 



(6.1) 



where Sq, p was defined by (4.2). 
De/me Woc/c diagonal R e M d 9mxdgm &y 



V 



(6.2) 



w/iere eac/i of R\, . . . , R m € Ri dx i d is itself block diagonal and given by 

V<Mn)- T 



Rt = ctt X 



(6.3) 



V0t(r,)" 



T/ie scalar a t used here was defined by (4.4). 

Finally, D £ ^gdmxgdm j s p 0S m ve definite diagonal and given by 



D 



\ 



D m J 



(6.4) 



where D t G M9 d>< 'J d 7 t — 1, . . . , m, is ^iven &y 



/ 0(0t(ri))det(V^(ri))wiJ 



V 



0(&(r a ))det(V&(r a ))w g J / 

(6.5) 

m which I denotes the d x d identity matrix. 
Remarks. Intuitively, this factorization of K decomposes finite element analysis into 
natural ingredients: A encodes the combinatorial connectivity of the mesh, R encodes 
the geometry of the mesh, S encodes the quadrature points, and D encodes the 
quadrature weights and conductivity. The scaling factor a t , which is necessary for 
our analysis, cancels out between D t and R t . 

Note that the positive definiteness of D follows from Assumptions 1, 2, and 3. 
Note also that A is a reduced node-arc incidence matrix of a multigraph defined on the 
nodes of T. (Multigraph indicates that more than one edge may connect a particular 
pair of vertices.) Each element gives rise to Z — 1 arcs in the graph. In particular, for 
each element t = 1, . . . , m, there is an arc joining each of its nodes 2, . . . , I to node 
1. Columns corresponding to nodes of Ti are omitted. (This is what is meant by 
"reduced.") 

Proof. Let £ be an arbitary vector in M™, and let u : f2 — > K be defined by 



1=1 
11 



For an element t, define function U t '■ T — » R by 

l 

Ut = J2^G M N^. (6.6) 

With these definitions, it is clear that u o (f> t = Ut- In (6.6), we follow the convention 
that & = in the case that i € {n + 1, . . . , n'}. 

The assembled stiffness matrix is defined by (3.7) so that 

n n 

i=l J=l 

n n m q 

= ^2^2&£ 1 ^2^2^xM ( f > t{rk)) ■ 6 l ((/)t(r fc ))V a; 7r :) (0 t (r fc ))det(V(/) t (r fc ))w fc 
»=i j=i t=i fe=i 

m q n n 

t=i k=i »=i j=i 

m <j 

= X! X! Va; ( rfc ) ) ■ ^ ( rfe ) ) v * ( rfc ) ) det ( ( rfc ) 

t=l fc=l 
m g 

= ^E]V x C/ t (r fe )^(^(r fe ))V x f/ t (r fe )det(V^(r fe )H 
t=i fe=i 

m 9 

= E E(V&(r fc )- r V,IMr fc )) • *(&(r fc )) 
t=i »=i 

(V0 t (r fc )- T V z (7 t (r fc ))det(V0 t (r fc )) Wfe (6.7) 
= v T L>v. (6.8) 

In (6.8), we have used the matrix D defined by (6.4). We have also introduced the 
vector v G M. dqm defined in block fashion as follows. Write v = [v\; ■ ■ ■ ;v m ] where 
v t £ R dg is itself composed of blocks v t = [i>t,i; • • • ; Here £ R d is defined to 
be 

vt.k = {V4>t{r k ))- T V z U t {r k )/a t . (6.9) 

It is clear by construction of D and v that v T Dv is equal to the expression in (6.7). 

Next, we claim that v = RSA£. Focus on the block corresponding to a particular 
element t £ {I, . . . , to} in which case we must show that vt = RtSQ tP A t £. The product 
A£ yields a vector with (I — I) entries that contains finite differences of entries of £. 
Specifically, the \i - 1 entry is &,G(t, M ) - £lg(m) for \i = 2, . . . , I. 

By definition of Sq, p in (4.2), it follows that the fc-block of entries (k = 1, . . . , q) 
of SQ,pA t £ is the d-vector \7 z U t (r k ), where 

Ut = (£.LG(t,2) - ^LG{t.l))N 2 H h {%LG(t,l) ~ ^LG(t,l)) N l- 

Comparing this equation to (6.6) indicates that U t and U t differ by 

^LG(t.l)(N 1 + --- + Ni). 

This latter quantity, however, is a constant function (because N\+- ■ -+Ni is identically 
1, a property that follows from (3.2) and the fact that I is a polynomial of degree at 
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most p and hence must be expressable as a sum of N^'s). Therefore, U t and U t have 
the same gradient. We conclude that the k block of Sq p A^ must equal V z (7 t (rfc). 
Combining this equation with (6.3) shows that 

RtS Q , p A t £, = (y4>t{r k )y T V z U t {r k )/a u 

and hence is equal to v t .k as defined by (6.9). This concludes the proof that v = RSA£. 

Since v T L>v = £ T K£, this means f ' A T S T R T D RS A£ = $, T K£. Since this holds 
for all £ £ R™ , and since a symmetric matrix C is uniquely determined by the mapping 
£ i ► £ T C£, this means that K = A T S T R T DRSA, which concludes the proof of the 
lemma. □ 

Now we analyze the singular values of J to finish the first part of Theorem 5.1. 
Let us define J = RS so that K = A T J T DJA as claimed. The block structures of R 
and S induce a corresponding block structure on J: 

/ Ji \ 

V Jm ) 

Jt = RtS Q . p (6.10) 



J 



where 



for all t = 1, . . . , to. Because of this structure, the maximum singular value of J is 
the maximum singular value among any of its blocks and similarly for its minimum 
singular value. Since in general cr max (AB) < cr max (y4)cr max (B), 

fmax(Jf) < 0'max(-Rt)Cmax(<S'Q,p) 

= cr m ax(diag(V0 t (ri)" T , . . . , V4>t(r q y T ))a QyP /a t 
= max cr max (V(/) t (r fc )~ 1 )cr Ci p/Q! t 

k=l,...,q 

= (TQ,p. (6.11) 

The last line follows from (4.4) and establishes (5.1). 

Since a m - ln (AB) > a m i n (A)a ln i n (B) for two matrices A, B with full column rank, 

= cr m in(diag(V0 t (ri)~ T , . . . , V0 t (rg)~ T ))T SiP /at 
= min CT m i n (V(/) t (r fc ) _1 )TQ.p/a t 

k—l,...,q 

= min (l/a max (V0t(rfc)))Tg )P /at 

k—l,...,q 

= • T Q, P / a t 

> t q , p / Ki (T). (6.12) 

For the last line, we used the fact that atfit < Ki(7"), which follows from (4.6). This 
establishes (5.2) and concludes the proof of the first factorization in Theorem 5.1. 

Our factorization K = A T J T DJ A is reminiscent of one proposed by Argyris [1] 
of the form K = ApA t , which he calls the "natural factorization." In Argyris's 
factorization, however, the matrix A has all +1 and entries and therefore is not 
a node-arc incidence matrix. The purpose of Argyris's matrix A is to assemble the 
element stiffness matrices, which constitute the blocks of the block-diagonal matrix 
P. 
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7. The second factorization. In this section we prove the second part of The- 
orem 5.1. We state the existence of the second factorization in the form of a lemma. 

Lemma 7.1. Let J and D be defined as in Lemma 6.1. Then J T DJ can be 
refactored as D 1 / 2 J T JD 1 / 2 , where 

• Matrix D e jm(l-i)xm(l-i) j s p 0S m ve definite diagonal and is written in 
block form 



D = m 



Q 



(7.1) 



where I is the (I — 1) x (I — 1) identity matrix. In this formula, ft is the 
minimum value of 9 over Gauss points of element t: 

ft = mm 9(Mrk)), (7.2) 

K=l,...,q 

and similarly, g t is the minimum value of det(V^t): 

g t = mm det(V^ t (r fe )). (7.3) 

k=l,...,q 

• Matrix J e R d 9mx(;-i)m is defined by 

J = D 1/2 JD~ 1/2 . (7.4) 

Proof. The fact that D 1 / 2 J T JD 1 / 2 = J T DJ follows as an immediate consequence 
of (7.4) regardless of how we have defined D. □ 

Next we analyze the condition number of H to finish proving the second part of 
Theorem 5.1. We define H = J T J so that J T DJ = B x l 2 Hb x l 2 . We must estimate 
the singular values of J, which are the square roots of the eigenvalues of H. 

Let the diagonal block of J associated with element t be denoted Jt for t = 
1, . . . , m. Because of the block structure, the maximum and minimum singular values 
for J are the maximum and minimum singular values among the blocks Jt which may 
be written 

1 rd/2 t j— 1/2 -1/2 -1/2 -l 

Jt = D t ' Jtft 9t m Q <*t 
= D t J t 

where 

A i— 1/2 -1/2 -1/2 -1 n l/2 
A = ft 9t m Q a t D t ■ 

Examining the constituent parts of D t given by (6.5), observing that the a 2 in (6.5) 
is cancelled by the a^ 1 in the preceding equation and applying the inequalities 

1 < LUk/m,Q < Mq/tuq 

for k = 1, . . . , q (see (4.1)), 



i < 0(<MnO)/t _1 < 



max fc=1 ,...^6>((/> t (r fc )) 
min fc= i ) ... i9 0(^(r fe )) 
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(see (7.2)), and 



i s a +fv7j. i \\ -i ^ max fe=i gdet(V0 t (r fc )) 

1 < det (V0 t (r fe ) 5t < — , , m t , 

min fe= i i ... i9 det(V0t(rfe)) 

(see (7.3)), we conclude that the diagonal entries of D t satisfy 

1 , a ,. .x . / max fc 6'(0 t (r fe ))\ 1/2 / max fc det(V<?!) t (r fe )) \ 1/2 / M s \ 1/2 
1 < D t (i,i) < ■ 



V min fe 6»(</> t (r fe ))/ \ min fc det(V^ t (r fc )) / \m 

By (4.7) and (4.8), the first two quantities on the right-hand side of the preceding 
equation are bounded by ^(T) 1 / 2 and ^(T) 1 / 2 respectively. Thus, it is apparent that 
for each i, 

1 < AM < (fl(T)«2(T)M e /m e ) 1/2 . 

Since J t = DtJt, we can combine the inequalities in the previous line with (6.11) and 
(6.12) to obtain: 

<W(J*) < HT^^T^M^aQjm 1 ^ 2 , (7.5) 

and 

Cmin(Jt) > tq iP /k\(T). (7.6) 

Since H = J T J and J has full column rank, (7.5) and (7.6) prove (5.3), which 
concludes the proof of Theorem 5.1. 

8. Preconditioning Strategy Summary. The main result of this paper is 
that the stiffness matrix K of (3.1) can be approximated by a symmetric diagonally 
dominant matrix K. In this section we discuss several approaches for using this 
approximation to efficiently solve the linear system K x = f by iteration. 

The first step in using the approximation to construct K. Recall that K — A T DA, 
and thus A and D must be formed. Forming A means construction of the multigraph 
consisting of I — 1 edges (a star-tree) per element of T. Matrix A is then the reduced 
node- arc incidence matrix of this multigraph as specified by Lemma 6.1. Construction 
of D is given by (7.1). Notice that computing A and D requires information about 
the mesh and original boundary value problem. In other words, our method is not 
applicable (at least not in a straightforward manner) if the only information about 
the original problem is the stiffness matrix K. 

Once K is on hand, there are several ways to proceed. The most straightforward 
is to suppose that one has an efficient preconditioned conjugate gradient solver for 
systems of the form if x = f . Efficient preconditioners for symmetric diagonally dom- 
inant systems were proposed and analyzed in a graph-theoretic framework by Vaidya. 
Vaidya's work is described and extended by [4, 6], and [11] contains a related analysis. 
One of the most recent improvements is due to Spiclman and Tcng [20] , who propose 
a graph-based preconditioner whose running time is 0(n 5 ' 4 ). (Spiclman and Tcng 
improve this bound to 0(n 1+£ ) for any e > 0, but the algorithm corresponding to 
the improved bound is no longer preconditioned conjugate gradients.) Other tech- 
niques proposed in the literature for symmetric diagonally dominant matrices such as 
algebraic multigrid could also be used. 
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If one has a good preconditioncr M for K, then one could also use M as a 
preconditioner directly for K. This is because of the "triangle inequality" (see [5]), 
which states 

a(K, M) < a(K, K)a(K, M). (8.1) 

This shows that the the overall support number is bounded by the product of the 
support numbers in each step of the approximation K w K w M . In particular. 
k(M, K) < k(H)k(M, K). Since the number of iterations of preconditioned conjugate 
gradients is bounded by the square root of the condition number, our analysis shows 
that the increase in the number of iterations for solving Kx = f (compared to solving 
Kx = f) is at most a factor of ^k(H), which in turn is bounded by 

e(T) 1/2 Ki(T)K 2 (T) 1/2 M e a ej ,/(maT ej ,). 

The asymptotically fastest known iterative algorithm for solving symmetric diag- 
onally dominant matrices is due to Elkin et al. [9]. This algorithm, which we denote 
EEST, extends Spiclman and Teng [20] and requires 0(n e (logn) s ) time to solve any 
diagonally dominant Kx — f , where n e is the number of nonzero entries in K and s 
is some constant. Our preconditioncr K is the Laplacian of a graph with ml edges, 
where m is the number of elements and I is the number of nodes in the reference 
element (refer to Lemma 6.1), hence n e = ml. We may assume that m = 0(n) for the 
following reasons. First, this relationship holds for all automatic finite element mesh 
generators of which we are aware. Second, a violation of the relationship m — 0(n) 
implies that the mesh has elements of unbounded aspect ratio which in turn means 
Ki(T) tends to infinity so that our results are not as useful anyway. Thus, it is safe to 
assume n e < const • nl, and we are regarding I (which depends on p, the polynomial 
degree, and d, the space dimension) as a constant. 

The EEST iteration is more complicated than preconditioned conjugate gradients, 
and hence it is unclear whether it is possible to simply replace the diagonally dominant 
coefficient matrix K by our stiffness matrix K and expect the algorithm to converge 
in 0(n(\ogn) s k(H) 1 / 2 ) time. 

On the other hand, one could obtain this running time 0{n{\ogn) s K{Hyl 2 ) by 
using a nested iteration: the outer iteration is preconditioned conjugate gradients 
in which K preconditions K. Thus, O^i/) 1 / 2 ) outer iterations are required. The 
inner loop is the EEST algorithm to apply the preconditioner, i.e., to solve Kx = f 
in 0(n(logn) s ) iterations. The EEST algorithm is fairly complex, and using it in 
a two-level manner would raise a number of difficulties associated with termination 
tests, so it is unlikely to be practical currently. 

This bound of 0(n(log n) s n(H) 1 ^ 2 ) on the number of operations is independent of 
the condition number of the underlying system because, as noted above k(H) depends 
on factors associated with the finite element method and on mesh quality measures 
but not on the conditioning of K . 

Our approximation scheme can be rewritten on an clcmcnt-by-element basis as 
follows. For some element index t G {1, . . . , m}, consider its element stiffness matrix 
K t — Aj Jj D t JfAt, where A t is defined in Lemma 6.1, Jt is defined by (6.10), and 
D t is defined by (6.5). The above proof shows that D t , which is the tth block of 
(7.1), is a good approximation to Jf D t Jt- Thus, we let K t = AjD t A t . The overall 
approximation is K = Y^itLi Ki- Note that the so-called splitting lemma frequently 
used in support-graph theory (see, e.g., [5]) shows that 

a(K, K) < maxiaiKuKi), a(K m , K m )), 
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Table 10.1 

Quadrature rules used for the six test cases. The rules are written in the form 
(ri, u>i), (r q ,w q ), where r^ is a Gauss point in the reference simplex and is its weight. In 
this table £1 = (10 - v / 20)/40 and £2 = 1 - 3£i. 



Test 


Q 


Ti 


((1/3, 1/3), 1/2) 


T 2 


((1/6, 1/6), 1/6), ((1/6, 2/3), 1/6), ((2/3, 1/6), 1/6) 


T 3 ,T 5 


((1/4, 1/4, 1/4), 1/6) 


T4, T 6 


((Ci, £1, €1), 1/24), ((&, £1, 6), 1/24), ((a, 6, £1), 1/24), (fe, £1, £1), 1/24) 



so it suffices to analyze the quality of approximation of -K"t to i^T t . 

Alternatively, we can take a global view as in the above analysis and write K — 
A T DA, where A = (Ai; A 2 ; . . . ; A m ) and D = diag(£»i, D 2 , D m ). To simplify 
notation, we have adopted the global view in this paper, but the reader should keep 
in mind that our approximation can take place element by element, which may be 
important in an implementation. 

9. Notes added in revision. Since the first version of this paper, there have 
been three related developments of interest. First, R. Gupta's Master's thesis [12] 
gave a geometric method to approximate element stiffness matrices by diagonally 
dominant matrices in the p = 1, d = 2 case. In follow up work, Wang and Sarin [23] 
showed how to parallelize the algorithm. Computational experiments are presented. 
For this particular case p = 1, d = 2, their result is essentially equivalent to ours. 

Second, Phillips and Miller [18] have shown that in the p = 1, d = 2 case, k(H) 
is still small even if some elements have very small angles, provided that no element 
has large angles (close to n). In our preceding analysis, the right-hand side of (5.3) 
would grow large in the presence of small angles since the factor K\(T) would be large. 
Wang and Sarin make a similar observation. 

In a third development, Avron ct al. [2] have shown how to get an approximation 
to element stiffness matrices that is optimal up to a constant factor. In other words, 
for an element stiffness matrix K t they find a diagonally dominant matrix K t such 
K(K t ,Kt) < c\K(K t ,Kl), where K' t is any other symmetric diagonally dominant ma- 
trix of the correct size and c\ depends only on p, d. Thus, their preconditioner could 
lead to a faster algorithm than ours since ours is not optimal in this sense. Their 
theory, however, does not subsume ours since they do not obtain any new bounds on 
k(K,K) comparable to (5.3). Their results are competitive and in some cases better 
than an algebraic multgrid method. 

10. Computational tests. There are at least two possible ways to test the 
effectiveness of our result: calculate the quality of matrix approximation or measure 
the speed of an iterative method. Since our theory is primarily about the former issue, 
and since Avron et al. conduct extensive testing on the latter question, we focus on 
the first question. 

We try out four specific commonly occurring cases: (d,p) = (2, 1), (2, 2), (3, 1), (3, 2),| 
i.e., linear and quadratic triangles and linear and quadratic tetrahedra. Six test 
meshes were tried, which we denote Ti, . . . , T$; the first two are two-dimensional and 
the last four are three-dimensional. For quadrature in the p = 1 cases, we use the 
midpoint rule. For the p = 2 cases, we use the symmetric d + 1-point rule, which is 
accurate for polynomials up to degree 2. Table 10.1 gives the quadrature rules used. 

Table 10.2 gives the values of relevant scalars for these test cases. It should be 
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Table 10.2 
Scalars relevant for test cases 



Test d p q o-Q^p tq iP MQ/m q 

Ti 2 I 1 LOO LOO LO 

T 2 2 2 3 5.26 0.83 1.0 

T 3 ,T 5 3 1 1 1.00 1.00 1.0 

T 4 ,T 6 3 2 4 6.47 0.63 1.0 



noted that <tq iP = tq, p = 1 for the p = 1 cases under the midpoint rule because in 
these cases Sq, p given by (4.2) turns out to be the identity matrix. 

Finally, we can tabulate element stiffness approximation bounds. (As mentioned 
in Section 8, it suffices to measure element stiffness approximation, which is much 
easier to compute than global stiffness approximation.) There are at least three rele- 
vant quantities to tabulate: xi — K {Kt,K t ), X2 = k(H), and X3 = 9(T)ki(T) 2 k 2 (T) • 
Maa f' p . which is the right-hand side of (5.3). Here, Kt is the stiffness matrix of 

element t (t = 1, . . . , m), and K t is the diagonally dominant preconditioner given by 
AjD t A t . It follows from Theorem 2.3 and (5.3) that xi < X2 < X3- 

We generate six meshes and take the max of Xi,X2,X3 over all elements of the 
mesh. For conductivity, we use 9 = 1, so that the factor 9(T) does not enter the 
bound. The meshes are generated as follows. In two dimensions, we use the Triangle 
mesh generator [19] to generate a mesh of an annulus, which is Ti and consists of 
linear (p = 1) elements. To generate T 2 , a quadratic mesh, we insert midpoint nodes 
in every element of Ti, and for the edges on the boundaries, we moved the midpoint 
nodes onto the boundary. 

The remaining four meshes are three-dimensional. For T3, we used the QMG 
mesh generator [17] to generate a mesh of a unit ball in R 3 . QMG's meshes are 
intended for the p = 1 case. For T4 we use midpoint insertion on the T3 mesh to 
obtain a p = 2 mesh. Again, midpoint nodes of edges whose endpoints were on the 
boundary were displaced outward onto the boundary. For T 5 and T 6 , we again mesh 
a unit ball (actually, one octant of the ball) using a mesh generator tailored for that 
domain only. The mesh generator for T 5 first generates a highly regular mesh of a 
unit tetrahedron, and then it projects the mesh nodes radially outward toward one 
facet so that an octant of the unit ball is covered. Finally, T 6 is obtained from T 5 
using midpoint insertion and displacement at the boundary. 

Then for every element in each mesh, we compute the three quantities Xi>X2>X3- 
We have tabulated the maximum values in Table 10.3. The undesirably large values 
of Xi for T3 and T4 appears to be due primarily to poorly shaped elements, i.e., large 
value of K\(T). This is evident from comparing T 3 and T4 against T 5 and T 6 , two 
meshes for the same domain but with much better shaped tetrahedra. 

The value of K2{T) was quite small in all tests; as mentioned earlier, it is iden- 
tically 1 for the p = 1 cases. For the p = 2 cases, interior elements are still linear, 
and elements adjacent to the boundary are fairly close to linear. It is interesting to 
note that xi = X2 in ah cases. This implies that the invariant subspaces of H corre- 
sponding to its extremal eigenvalues meet the range space of D 1 / 2 A. We do not have 
a deeper explanation for this observation. 

Although an iterative method was not tested, some conclusions can still be drawn 
from Table 10.3. For example, if our approximation were used on Tq in an iterative 
setting, then the slowdown would be at most a factor of \J K(K t ,K t ), i.e., a slowdown 
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Table 10.3 

Approximation of the element stiffness matrix by a diagonally dominant matrix. 



Test 


m 


n 




*2(T) 


Xi 


X2 


X3 


Ti 


234 


143 


4.1 


1.0 


16.6 


16.6 


16.6 


T 2 


234 


520 


4.7 


1.3 


185.1 


185.1 


1117.2 


T 3 


9078 


1913 


59.6 


1.0 


3549.5 


3549.5 


3549.5 


T 4 


9078 


13608 


90.2 


2.7 


5.96- 10 4 


5.96- 10 4 


2.32 • 10 6 


T 5 


729 


220 


5.0 


1.0 


24.9 


24.9 


24.9 


T 6 


729 


1330 


5.0 


1.2 


396.6 


396.6 


2644.5 



of at most a factor of 20. In practice, the method seems to work better than that 
according to Avron et al. 

11. Open questions. This work is the first to extend support-tree methods, 
which previously have been shown to be good preconditioncrs for diagonally dominant 
matrices with negative off-diagonal entries, to the class of finite clement matrices. We 
have shown that the scope of the method includes a standard scalar elliptic boundary 
value problem, but perhaps the scope of finite element problems that can be tackled 
with this method could be expanded further. 

One generalization would be the class of problems V-(0(a;)Vti) = — /, where Q(x) 
is a spatially varying d x d symmetric positive definite matrix. This generalization 
would present problems for our current analysis in the case that O(x) is highly ill- 
conditioned. It would still be straightforward to write K — A T J T DJA where D is 
now block diagonal, but our analysis of the introduction of J would run into trouble 
because the dq x dq diagonal blocks of D are no longer individually well conditioned. 

It would also be interesting to tackle vector problems such as linear elasticity or 
Stokes' flow, or higher-order equations like the biharmonic equation. It seems likely 
that our techniques can extend to at least some of these problems since they all have a 
symmetric positive definite weak form. On the other hand, a recent result [7] suggests 
that the extension of our results to linear elasticity will not be straightforward because 
of the higher nullity of element stiffness matrices in the case of linear elasticity. A fur- 
ther generalization would be to unsymmetric problems like the convection-diffusion 
equation. The latter class of problems would require substantial rethinking of the 
whole approach since condition number reduction, which is very relevant for the ap- 
plication of conjugate gradients to symmetric positive definite systems, is less relevant 
to the application of GMRES to unsymmetric systems. 

Our analysis is based on condition numbers (support numbers). One drawback 
of this approach is that the convergence and work estimates may be too pessimistic. 
For instance, the condition number of the preconditioned linear systems depends 
on Ki(T), the worst aspect ratio of any element in the mesh. If there is only one 
poorly shaped element in the mesh, we expect iterative solvers will only take a few 
extra iterations since changing a single element implies a low-rank correction to the 
assembled stiffness matrix. Any analysis based on condition numbers will be unable 
to capture this effect. A related open issue is whether we can exploit recent work in 
mesh quality metrics [16] to show that "good meshes" both have small error in the 
FEM approximation and also produce linear systems that can be well approximated 
by diagonally dominant systems. 

Another point to make about our method is that, although the condition number 
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of the preconditioned system has an upper bound independent of 

Rq = max #(2;)/ min#(a;), 

there will still be a loss of significant digits due to roundoff error when using our 
method in the case that Rg is large. This is because the system matrix and the 
preconditioner separately are ill-conditioned operators. A special case of (3.1) in 
which Rg is extremely large and in which / = was considered in [22]. That paper 
proposed a method based on Gaussian elimination in which the loss of significant 
digits is avoided. Some of the ideas behind [22] were also extended to solution via 
conjugate gradient using support preconditioners by [14]. The methodology in [14], 
however, was for node-arc incidence matrices and for a particular kind of support 
preconditioner called a support tree [11], and it is not clear whether that method 
would apply to the present setting. 

Finally, there is greater effort to couple automatic finite clement mesh generation 
to detailed analysis of the solver and its error bounds. Our approximation result adds 
another ingredient to this mix; it would be interesting if an automatic finite element 
mesh generator could be tailored to reduce specifically the quantity k(K, K). 
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